The author acknowledges the assistance of AI tools in the preparation of this document (GitHub Inc. 2025; OpenAI 2025).
Loading Data
Code
# Libraries
library(readxl)
library(dplyr)
library(knitr)
library(kableExtra)
library(moments)
library(ggplot2)
suppressWarnings(library(plotly))
library(lme4)
library(performance)
library(sjPlot)
library(lmerTest)
library(DT)
# Excel data originating from the biodiversity_abundance data sheet
mydata <- read_excel("Data.xlsx")
# Convert categorical variables to factors
mydata$Month <- factor(mydata$Month)
mydata$Hive <- factor(mydata$Hive)Code
# Data table
datatable(
mydata,
options = list(scrollX = TRUE),
rownames = FALSE
)Exploratory Data Analysis
Code
# Summary statistics for chill coma onset time by hive and month
summary <- mydata %>%
group_by(Hive, Month) %>%
summarise(
Mean = mean(`Coma time`),
Median = median(`Coma time`),
IQR = IQR(`Coma time`),
SD = sd(`Coma time`),
SE = SD / sqrt(n()), # Standard error
Skewness = skewness(`Coma time`),
Kurtosis = kurtosis(`Coma time`),
n = n(),
.groups = "drop"
)
# Summary statistics for chill coma onset time by month
overall_summary <- mydata %>%
group_by(Month) %>%
summarise(
Hive = "Overall",
Mean = mean(`Coma time`),
Median = median(`Coma time`),
IQR = IQR(`Coma time`),
SD = sd(`Coma time`),
SE = SD / sqrt(n()), # Standard error
Skewness = skewness(`Coma time`),
Kurtosis = kurtosis(`Coma time`),
n = n(),
.groups = "drop"
)
# Combine and arrange
full_summary <- bind_rows(summary, overall_summary) %>%
arrange(Month, Hive)
# Render table
kable(full_summary, digits = 2)| Hive | Month | Mean | Median | IQR | SD | SE | Skewness | Kurtosis | n |
|---|---|---|---|---|---|---|---|---|---|
| Hive 1 | April | 107.20 | 113.37 | 27.97 | 21.13 | 6.68 | -0.07 | 2.07 | 10 |
| Hive 2 | April | 122.22 | 119.36 | 14.73 | 27.66 | 8.75 | 1.25 | 4.75 | 10 |
| Hive 3 | April | 108.52 | 111.88 | 20.63 | 19.00 | 6.01 | -1.27 | 4.17 | 10 |
| Overall | April | 112.65 | 114.30 | 25.55 | 23.15 | 4.23 | 0.69 | 5.61 | 30 |
| Hive 1 | May | 141.68 | 145.20 | 30.71 | 28.17 | 8.91 | -0.19 | 2.29 | 10 |
| Hive 2 | May | 148.28 | 152.70 | 26.00 | 27.05 | 8.56 | -1.50 | 4.77 | 10 |
| Hive 3 | May | 189.27 | 164.82 | 85.41 | 74.69 | 23.62 | 0.59 | 2.21 | 10 |
| Overall | May | 159.74 | 154.02 | 30.85 | 51.61 | 9.42 | 1.49 | 5.74 | 30 |
Code
colour_scheme <- c("April" = "darkgreen", "May" = "lightgreen")
summary_df <- mydata %>%
group_by(Month) %>%
summarise(
mean = mean(`Coma time`),
se = sd(`Coma time`) / sqrt(n())
)
# Column graph
graph <- ggplot(summary_df, aes(x = Month, y = mean, fill = Month)) +
geom_col(width = 0.4, color = "black") +
geom_errorbar(aes(ymin = mean - se, ymax = mean + se), width = 0.2) +
labs(
y = "Chill Coma Onset Time (s)",
x = "Month"
) +
scale_fill_manual(values = colour_scheme) +
scale_y_continuous(expand = expansion(mult = c(0, 0.05)), limits = c(0, NA)) + # y-axis starts at zero
theme_minimal(base_size = 14) +
theme(
panel.grid = element_blank(),
axis.line.x = element_line(color = "black", linewidth = 0.7),
axis.line.y = element_line(color = "black", linewidth = 0.7),
axis.ticks = element_line(color = "black"),
legend.position = "right"
)
ggplotly(graph)Code
# Boxplot
boxplot <- mydata %>%
ggplot(aes(x = Hive, y = `Coma time`, fill = `Month`)) +
geom_boxplot(position = position_dodge(width = 0.8)) +
labs(
x = "Hive ID",
y = "Chill Coma Onset Time (s)",
) + scale_fill_manual(values = c("April" = "darkgreen", "May" = "lightgreen")) + theme(
panel.grid = element_blank(),
axis.line.x = element_line(color = "black", linewidth = 0.7),
axis.line.y = element_line(color = "black", linewidth = 0.7),
axis.ticks = element_line(color = "black"),
legend.position = "right")
ggplotly(boxplot) %>%
layout(
plot_bgcolor = "white", # plot area background
paper_bgcolor = "white" # overall background
)Code
# Histogram
month_colors <- c("April" = "darkgreen", "May" = "lightgreen")
histogram <- mydata %>%
ggplot(aes(x = `Coma time`, fill = Month)) +
geom_histogram(color = "black", bins = 10, position = "identity", alpha = 0.7) +
labs(x = "Chill Coma Onset Time (s)", y = "Frequency") +
facet_wrap(~ Hive, scales = "free_x") +
scale_fill_manual(values = month_colors) +
theme_minimal() + theme(
panel.grid = element_blank(),
axis.line.x = element_line(color = "black", linewidth = 0.7),
axis.line.y = element_line(color = "black", linewidth = 0.7),
axis.ticks = element_line(color = "black"),
legend.position = "right")
# Display the plot
histogramAssumptions
Code
# Model assumptions
model <- lmer(`Coma time` ~ Month + (1 | Hive), data = mydata)
performance::check_model(model)Code
# Model assumptions for square-root transformation
mydata$Coma_time_sqrt <- sqrt(mydata$`Coma time`)
model_sqrt <- lmer(Coma_time_sqrt ~ Month + (1 | Hive), data = mydata)
performance::check_model(model_sqrt)Model
Code
# Model
model_lmerTest <- lmer(`Coma time` ~ Month + (1 | Hive), data = mydata)
tab_model(
model_lmerTest,
show.se = TRUE,
show.stat = TRUE,
show.ci = FALSE,
show.re.var = TRUE,
dv.labels = "Chill Coma Onset time (s)",
string.pred = "Predictor",
string.est = "Estimate"
)| Chill Coma Onset time (s) | ||||
|---|---|---|---|---|
| Predictor | Estimate | std. Error | Statistic | p |
| (Intercept) | 112.65 | 8.71 | 12.93 | <0.001 |
| Month [May] | 47.09 | 10.16 | 4.63 | <0.001 |
| Random Effects | ||||
| σ2 | 1549.42 | |||
| τ00 Hive | 72.73 | |||
| ICC | 0.04 | |||
| N Hive | 3 | |||
| Observations | 60 | |||
| Marginal R2 / Conditional R2 | 0.258 / 0.291 | |||
References
GitHub Inc. 2025, GitHub Copilot [Code completion tool], GitHub Inc., https://github.com/copilot.
Native Bee Hives (n.d.) Tetragonula carbonaria [Photograph]. Available at: https://www.nativebeehives.com/wp-content/uploads/photo-gallery/nativebeedaisy1.jpg (Accessed: 27 May 2025).
OpenAI 2025, ChatGPT (Version 4), OpenAI, viewed 18 May 2025, https://openai.com/chatgpt.